#include "AnaDef.h"
#include "XGLUtils.h"
#include "AliTRDdEdxBaseUtils.h"

void drawBG(const TString ilist, const Int_t v0pid)
{
  TList * ll = new TList;

  TCanvas *c1=new TCanvas;
 
  TChain *ch=XGLUtils::InputFiles(ilist,"tree");
  const Int_t nch = 6;
  const Int_t ncls = 18;

  if(v0pid!=2 && v0pid!=4){
    printf("wrong v0pid! %d\n", v0pid);exit(1);
  }

  const TString qa = Form("tpcncls>120 && %s", v0pid==2? "(pidV0==1 || pidV0==2)" : "pidV0==4");

  const TString meancut = Form("trdncls/trdnch>=%d && trdnch>=%d && %s", ncls,  nch, qa.Data());
  const TString tag(Form("nch%d_ncls%d_pid%d",nch, ncls, v0pid));

  //===============================
  TH2D *h1=new TH2D(Form("hsigvsbg%d%d", nch, ncls), meancut, 100,0.3, 1e4, 200,0, 10); ll->Add(h1);
  AliTRDdEdxBaseUtils::BinLogX(h1->GetXaxis());

  h1->SetMinimum(1);

  //pidTPC0 does not help
  ch->Draw(Form("trdsig:bg>>%s", h1->GetName()), h1->GetTitle(), "colz");

  //=============================
  TH1D *hnor=0, *hmpv=0, *hwid=0, *hres=0, *hchi=0;

  TList * lslice = new TList;
  XGLUtils::FitSlicesY(h1, XGLUtils::RawGaus, hnor, hmpv, hwid, hres, hchi, 50,  -999, -999, lslice);

  ll->Add(hnor);ll->Add(hmpv);ll->Add(hwid);ll->Add(hres);ll->Add(hchi);
  //======================
  TH2D *h2=XGLUtils::NormalHist(h1, 50,1); ll->Add(h2);

  //======================
  Double_t pars[]={0.72, 1.33, 7.7, 0.313, 3.05, 0.035, 2.17, 0.66};
    
  Double_t errs[100], chi[100];
  XGLUtils::SetFitPrintLevel(1);
  const Int_t kfail = XGLUtils::ChisquareFit(hmpv, AliTRDdEdxBaseUtils::MeandEdxTR, 8, pars, errs, chi);
 
  TH1D *hfit = XGLUtils::GetHfit("hfit", AliTRDdEdxBaseUtils::MeandEdxTR, pars, 0.001, 4e4, kTRUE); ll->Add(hfit);

  TF1 * ffit = new TF1("ffit", AliTRDdEdxBaseUtils::MeandEdxTR, 0.1, 1e5, 8); ll->Add(ffit);
  ffit->SetParameters(pars);
  ffit->SetChisquare(chi[0]);
  ffit->SetNDF(chi[1]);
 
  if(kfail){
    hfit->SetTitle(Form("Fit fail %d!\n", kfail));
  }
  else{
    hfit->SetTitle(Form("%.3e , %.3e , %.3e , %.3e , %.3e , %.3e , %.3e , %.3e , (%.2f/%d) %s %s", ffit->GetParameter(0), ffit->GetParameter(1), ffit->GetParameter(2), ffit->GetParameter(3), ffit->GetParameter(4), ffit->GetParameter(5), ffit->GetParameter(6), ffit->GetParameter(7), ffit->GetChisquare(), ffit->GetNDF(), ilist.Data(), hmpv->GetTitle()));
  }

  //TF1::Fit is not stable Warning in <Fit>: Abnormal termination of minimization.
  /*
  TF1 * ffit = new TF1("ffit", AliTRDdEdxBaseUtils::MeandEdxTR, 0.1, 1e5, 8);
  ffit->SetParameters(pars);
  printf("fitting ffit...\n\n");
  hmpv->Fit(ffit);
  */

  //========================
  gPad->SetLogz(0); 
  gPad->SetLogx(1); 
    
  h2->Draw("colz"); 
    
  ffit->SetLineColor(kBlack);
  //ffit->SetLineWidth(2);
  ffit->Draw("same");
    
  c1->Print(tag+"_2d.eps");

  //===============
    
  hmpv->SetMarkerStyle(29);
  hmpv->SetLineColor(kBlack);
  hmpv->Draw();
    
  ffit->SetLineColor(kRed);
  //ffit->SetLineWidth(2);
  ffit->Draw("same");

  c1->Print(tag+"_1d.eps");

  //========================>>
  gPad->SetLogx(0);
  gPad->SetLogz(0);

  const TString widcut = qa;

  printf("widcut for scaled: %s\n", widcut.Data());

  //scaled signal
  TH2D * hscaled = new TH2D("hscaled", widcut, 100, 0.3, 1e4, 200, 0,2); ll->Add(hscaled);
  AliTRDdEdxBaseUtils::BinLogX(hscaled->GetXaxis());

  ch->Draw("trdsig/expmpv:bg>>hscaled", hscaled->GetTitle(), "colz");
  hscaled->Draw("colz");
  c1->Print(tag+"_scaledbg.eps");

  //--
  TH2D * hscn0 = new TH2D("hscn0", widcut,130, 10, 140, 200, 0,2); ll->Add(hscn0);

  ch->Draw("trdsig/expmpv:trdncls>>hscn0",hscn0->GetTitle(), "colz");
  hscn0->Draw("colz");
  c1->Print(tag+"_scaledN.eps");

  //--

  TH1D *hSanor=0, *hSampv=0, *hSawid=0, *hSares=0, *hSachi=0;
  //has to use func fit, otherwise there is not error for width and TF1::Fit will say "Fit data is empty"
  XGLUtils::FitSlicesY(hscn0,  XGLUtils::RawGaus, hSanor, hSampv, hSawid, hSares, hSachi, 50, -999, -999, lslice);

  ll->Add(hSanor);ll->Add(hSampv);ll->Add(hSawid);ll->Add(hSares);ll->Add(hSachi);
  //--

  TH1D *hfinal = hSawid; 
  hfinal->Draw("e");

  TF1 * finv = new TF1("finv","[0]+[1]*pow(x,[2])", 10, 150); ll->Add(finv);
  finv->SetLineColor(kRed);
  finv->SetLineWidth(1);
  finv->SetParameters(0, 1.451, -0.5);
  printf("fitting width...\n");
  hfinal->Fit(finv);

  hfinal->Draw("same");

  hfinal->SetTitle(Form("%.3e , %.3e , %.3e , (%.2f/%d) %s %s", finv->GetParameter(0), finv->GetParameter(1), finv->GetParameter(2), finv->GetChisquare(), finv->GetNDF(), ilist.Data(), hfinal->GetTitle()));

  TLegend *reslg = new TLegend(0.3, 0.7, 0.9, 0.9);
  reslg->AddEntry(hfinal,"p, #pi, e (p-Pb #sqrt{s_{NN}}=5.02 TeV)","p");
  reslg->AddEntry(finv,Form("fit, %s %.2f %.2f %.2f", finv->GetTitle(),  finv->GetParameter(0), finv->GetParameter(1), finv->GetParameter(2)));
  reslg->Draw();

  c1->Print(tag+"_resolution.eps");

  //========================<<

  TFile *fout=new TFile(Form("calib_%s.root", tag.Data()),"recreate");
  ll->Write();
  lslice->Write("slice", TObject::kSingleKey);
  fout->Save();
  fout->Close();
}


int main(int argc, char *argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=3){
    printf("argc!=2\n"); return 1;
  }

  drawBG(argv[1], atoi(argv[2]));
  return 0;
}
